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Abstract:  A  new  material  model  for  sand  has  been  developed  in  order  to  include  the  effects  of 
the  deformation  rate  and  the  degree  of  saturation  on  the  constitutive  response  of  this  material. 
The  model  is  an  extension  of  the  original  high  strain-rate  compaction  model  for  sand  developed 
by  Laine  and  Sandvik  and  an  elastic-visco-plastic  material  model  for  sand  recently  proposed 
by  Tong  and  Tuan  in  which  these  effects  were  neglected.  The  new  material  model  was  parame¬ 
terized  using  the  available  experimental  data  for  sand  with  different  levels  of  saturation  tested 
mechanically  at  different  strain  rates.  The  model  is  next  used,  within  a  non-linear-dynamics  tran¬ 
sient  computational  analysis,  to  study:  (a)  various  phenomena  associated  with  the  explosion  of 
shallow- buried  and  ground-laid  mines  and  (b)  the  dynamic  behaviour  of  a  vehicle  during  an  off¬ 
road  ride.  The  computational  results  are  then  compared  with  the  corresponding  experimental 
results.  This  comparison  suggested  that  the  newly  developed  material  model  for  sand  captures 
the  essential  features  of  the  dynamic  behaviour  of  sand  with  different  levels  of  saturation  when 
subjected  to  realistic  high  and  low  strain-rate  loading  conditions. 
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1  INTRODUCTION 

Recent  advances  in  numerical-analysis  capabilities, 
particularly  the  coupling  of  Eulerian  solvers  (used 
to  model  gaseous  detonation  products  and  air) 
and  Lagrangian  solvers  (used  to  represent  vehi¬ 
cles/platforms,  tires,  and  soil),  have  allowed  simula¬ 
tions  to  provide  insight  into  complex  loading  created 
by  the  detonation  of  buried  landmines  as  well  as  the 
assessment  of  the  off-road  (uneven-terrain)  vehicle 
dynamics  (including  vehicle  rollover  stability,  crew 
and  on-board  equipment  survivability,  etc.).  How¬ 
ever,  a  quantified  understanding,  through  computer 
modelling,  of  both  the  landmine-detonation-induced 
loading  and  the  resulting  interactions  between  the 
detonation  products,  mine  fragments  and  soil  ejecta 
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with  the  target  structures /vehicles  as  well  as  of  the 
off-road  vehicle  dynamics  is  still  not  mature.  As 
discussed  in  our  previous  work  [1-3],  the  lack  of 
maturity  of  computer  simulations  of  the  landmine- 
detonation  events  is  mainly  due  to  inability  of  the  cur¬ 
rently  available  material  models  to  realistically  repre¬ 
sent  the  response  of  the  materials  involved  (primarily 
soil)  under  high-deformation,  high-deformation-rate, 
and  high-temperature  conditions,  and  the  type  of  con¬ 
ditions  accompanying  a  landmine  detonation.  Similar 
deficiencies  (primarily  those  related  to  a  lack  of  con¬ 
sideration  of  the  effect  of  moisture  and  various  organic 
and  inorganic  constituents)  in  the  soil-material  mod¬ 
els,  but  under  lower- deformation  rates,  are,  at  least 
partly,  responsible  for  inadequacies  of  the  computer- 
based  off-road  dynamics  simulations  [4].  This  is  the 
main  reason  that  the  objective  of  the  present  work  was 
to  attempt  to  develop  a  visco-plastic  material  model 
for  soil  that  provides  a  proper  account  for  the  physi¬ 
cal  and  chemical  make-up  of  the  soil  and  is  suitable 
for  both  mine-blast  and  off-road  vehicle  dynamics 
computational  analyses. 
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In  the  present  work,  the  problem  of  material  model 
derivation  and  validation  for  ‘cohesion  less’  sand- 
based  soils  at  various  saturation  levels  is  addressed. 
Since  this  model  was  jointly  developed  by  Clemson 
University  (CU)  and  the  Army  Research  Laboratory 
(ARL),  Aberdeen,  Proving  Ground,  MD,  it  will  be 
referred  to,  in  the  remainder  of  the  document,  as  the 
CU-ARL  visco-plastic  sand  model. 

Soil  is  a  very  complex  material  whose  properties 
vary  greatly  with  the  presence /absence  and  relative 
amounts  of  various  constituent  materials  (sand,  clay, 
silt,  gravel,  etc.),  particle  sizes,  and  particle  size  distri¬ 
bution  of  the  constituent  materials.  In  addition,  the 
moisture  content  and  the  extent  of  precompaction 
can  profoundly  affect  the  soil  properties.  In  the  mine- 
blast  computational- analysis  community,  the  so- 
called  ‘porous-material/compaction’  model  proposed 
by  Laine  and  Sandvik  [5]  has  been,  for  quite  some 
time,  the  soil  model  that  provided  the  best  compro¬ 
mise  between  the  inclusion  of  essential  physical  phe¬ 
nomena  reflecting  material  response  under  dynamic 
loading  and  computational  simplicity.  However,  this 
model  was  developed  essentially  for  dry  sand  and,  as 
demonstrated  by  many  researchers  (e.g.  references  [6] 
to  [8]),  cannot  account  for  the  effects  of  moisture,  clay, 
and/or  gravel  in  soil.  To  overcome  these  deficiencies 
of  the  original  porous-material/compaction  model, 
CU  and  the  ARL  jointly  developed  [9-13]  and  subse¬ 
quently  parameterized  (using  the  results  of  a  detailed 
investigation  of  dynamic  response  of  sand  at  differ¬ 
ent  saturation  levels,  as  carried  out  by  researchers  at 
the  Cavendish  Laboratory,  Cambridge,  UK  [14,  15])  a 
series  of  CU-ARL  sand  models  capable  of  capturing 
the  effect  of  moisture,  clay  and  gravel. 

The  aforementioned  CU-ARL  soil  material  mod¬ 
els  were  all  rate-independent  elastic-plastic  in  nature 
and  optimized  with  respect  to  the  large-strain/high- 
deformation-rate/high-pressure  conditions  encoun¬ 
tered  during  detonation  of  landmines  buried  in  soil. 
As  such  these  models  are  not  suitable  for  use  in  the 
analysis  of  off-road  tire-soil  interactions  (critical  for 
the  rough-terrain  vehicle  dynamics  and  stability  anal¬ 
yses).  Consequently,  an  attempt  is  made  in  the  present 
work  to  use  the  basic  understanding  of  sand  dynamic 
behaviour  and  data  acquired  during  the  development 
of  the  CU-ARL  sand  model  in  order  to  extend  this 
model  to  the  deformation  regimes  encountered  during 
off-road  tire-soil  interactions. 

Sand  has  generally  a  complex  structure  consisting  of 
mineral  solid  particles  that  form  a  skeleton.  The  rela¬ 
tive  volume  fractions  of  the  three  constituent  materials 
in  the  sand  (the  solid  mineral  particles,  water,  and  air) 
are  generally  quantified  by  the  porosity,  apor,  and  the 
degree  of  saturation  (saturation  ratio),  /3sat,  which  are 
respectively  defined  as 

-  Y?  m 

^por  —  x  r  (1) 


and 

ftat  =  ^  (2) 

vp 

where  Vp  is  the  volume  of  voids  (pores),  Vw  is  the 
volume  of  water,  and  V  is  the  total  volume. 

Deformation  of  the  sand  is  generally  believed  to 
involve  two  main  basic  mechanisms  [16, 17]:  (a)  elastic 
deformation  (at  low-pressure  levels)  and  fracture  (at 
high-pressure  levels)  of  the  inter-particle  bonds  and 
(b)  elastic  and  plastic  deformations  of  the  three  con¬ 
stituent  materials  (sand  particle,  water,  and  air)  in  the 
sand.  The  relative  contributions  of  these  two  defor¬ 
mation  mechanisms  as  well  as  their  behaviour  are 
affected  primarily  by  the  degree  of  saturation  of  sand 
and  the  deformation  rate.  Specifically,  in  dry  sand, 
the  first  mechanism  controls  the  sand  deformation  at 
low  pressures  while  the  second  mechanism  is  domi¬ 
nant  at  high  pressures  and  the  effect  of  deformation 
rate  is  of  a  second  order.  In  sharp  contrast,  in  satu¬ 
rated  sand,  very  low  inter- particle  friction  diminishes 
the  role  of  the  first  deformation  mechanism.  On  the 
other  hand,  the  rate  of  deformation  plays  an  important 
role.  At  low  deformation  rates,  the  water/ air  residing 
in  the  sand  pores  is  squeezed  out  during  deforma¬ 
tion  and,  consequently,  the  deformation  of  the  sand  is 
controlled  by  the  deformation  of  the  solid  mineral  par¬ 
ticles.  At  high  pressures,  on  the  other  hand,  water/ air 
is  trapped  within  the  sand  pores  and  the  deforma¬ 
tion  of  the  sand  is  controlled  by  the  deformation  and 
the  volume  fractions  of  each  of  the  three  constituent 
phases.  Thus,  for  a  material  model  for  sand,  such  as 
the  CU-ARL  visco-plastic  sand  model  developed  in  the 
present  work,  to  be  considered  as  realistic  (i.e.  of  high 
fidelity),  it  should  be  able  to  reflect  these  differences  in 
the  mechanical  response  of  sand  under  different  pres¬ 
sure  and  strain-rate  loading  conditions.  In  addition, 
removal  of  the  main  deficiencies  of  the  so-called  ‘effec¬ 
tive  stress  approach’  [16-21]  (i.e.  the  inability  of  this 
approach  to  account  for:  (a)  deformation  of  the  solid 
particles  under  shock  loads  and  (b)  for  the  fact  that 
due  to  a  very  short  duration  of  shock  loading,  water 
may  become  trapped  in  soil  pores  and  provide  addi¬ 
tional  load  support)  is  also  emphasized  in  the  present 
work. 

The  organization  of  the  article  is  as  follows.  A 
brief  overview  of  the  prior  elastic-plastic  and  elastic- 
visco-plastic  material  models  for  sand  is  provided  in 
section  2.1.  Determination  of  the  material  evolution 
during  loading  and  numerical  solution  of  the  associ¬ 
ated  governing  equations  are  discussed  in  sections  2.2 
and  2.3,  respectively.  Derivation  and  parameterization 
of  the  newly  developed  CU-ARL  visco-plastic  sand 
material  model  are  presented  in  section  2.4.  Validation 
of  the  same  model  against  a  set  of  high-deformation- 
rate  (landmine  blast)  and  low- deformation  rate  (off¬ 
road  vehicle  dynamics)  results  is  given  in  section  3. 
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A  brief  summary  and  the  conclusions  obtained  in  the 
present  work  are  discussed  in  section  4. 

2  MODEL  DEVELOPMENT  AND 
COMPUTATIONAL  PROCEDURE 

In  the  following  subsections  of  the  present  section, 
a  description  will  be  provided  of  some  earlier  sand 
material  models  and  of  the  CU-ARL  visco-plastic 
sand  model.  It  should  be  first  recognized  that  since 
these  models  are  of  an  either  elastic-plastic  or  elastic- 
visco-plastic  nature,  they  are  composed  of  the  follow¬ 
ing  three  main  components:  (a)  a  yield  function  that 
defines  the  stress  based  condition(s)  that  must  be  sat¬ 
isfied  before  plastic-visco-plastic  yielding  can  begin; 
(b)  a  material  constitutive  relation  that  governs  the 
evolution  of  the  yield  function  during  elastic-visco- 
plastic  deformation;  and  (c)  a  flow  rule  describing 
the  evolution  of  the  plastic-strain  components  during 
plastic  deformation. 

2. 1  Prior  sand  modelling  efforts 

2.1.1  Yield  criterion 

In  most  early  investigations  of  mine  blast  and  off-road 
vehicle  dynamics  problems,  the  mechanical  response 
of  soils  was  modelled  using  the  rate-independent 
elastic-ideal-plastic  Drucker-Prager  model  [21,  22] 
within  which  the  yield  function  is  defined  as 

/  (J21  h)  —  y/Ji  —  a  +  0I\  (3) 

where  and  J2  are  the  first  and  the  second  invariant 
of  the  total  stress  and  the  stress  deviator,  respectively, 
while  a  is  the  zero-pressure  (shear)  cohesive  strength 
parameter  and  0  is  a  friction-angle-dependent  coef¬ 
ficient.  Plastic  yielding  occurs  when  f(j2,h)^0. 
It  should  be  noted  that,  equation  (3)  reduces  to 
the  standard  von  Mises  yield  function  when  0  =  0. 
Also,  equation  (3)  simply  accounts  for  the  fact  that 
inter-particle  cohesion  imparts  a  zero-pressure  shear 
strength  (a)  to  the  soil,  while  the  application  of  con¬ 
fining  pressure  (P  =  — fi/3)  causes  the  effective  shear 
strength  of  soil  to  increase  (due  to  an  increase  in 
the  inter-particle  friction) .  It  should  be  further  noted 
that,  throughout  this  article,  the  following  stress  sign 
convention  is  used:  tensile  stresses  are  assigned  a  pos¬ 
itive  value  while  compressive  stresses  are  negative.  A 
schematic  of  the  Drucker-Prager  model,  equation  (3), 
is  given  in  Fig.  1  (a) . 

Early  applications  of  the  original  Drucker-Prager 
model  (e.g.  references  [22]  and  [23])  in  various  soil- 
based  structural  analyses  revealed  some  serious  short¬ 
comings  of  this  model. 

1.  The  extent  of  plastic  dilatancy  (i.e.  the  increase 
in  volume)  under  shear  loading  is  typically  much 


(b) 

Fig.  1  Simple  schematic  of:  (a)  the  original  Drucker- 
Prager  elastic-ideal-plastic  yield  function  and 
(b)  the  modified  Drucker-Prager  model  with  a 
plastic -volume  hardening-based  moving  cap 

greater  than  that  observed  experimentally  (in  fact 
in  some  cases,  shear  loading  gives  rise  to  a  decrease 
in  volume). 

2.  Hydrostatic  loading/unloading  is  generally  found 
to  be  accompanied  by  a  considerable  hysteresis  that 
cannot  be  accounted  for  by  the  original  Drucker- 
Prager  model  equation  (3) ,  in  which  yield  function 
does  not  intersect  the  hydrostatic-loading  fi  axis 
(at  high  pressures)  and  when  the  hydrostatic  elastic 
response  is  associated  with  a  constant  value  of  the 
bulk  modulus. 

3.  At  very  high  pressure  levels,  soil  tends  to  behave 
more  as  a  fluid-like  material.  Under  such  condi¬ 
tions,  the  shear  strength  of  the  material  is  not 
expected  to  vary  with  hydrostatic  pressure.  In  other 
words,  the  yield  function  should  be  essentially  inde¬ 
pendent  of  7i  at  high  pressures.  According  to  the 
Drucker-Prager  yield  function,  equation  (3),  how¬ 
ever,  shear  strength  of  the  soil  continues  to  increase 
indefinitely  with  an  increase  in  pressure. 

To  overcome  the  aforementioned  shortcomings 
of  the  Drucker-Prager  elastic-ideal-plastic  model, 
Dimaggio  and  Sandler  [23],  Sandler  et  al.  [24], 
and  Sandler  and  Rubin  [25]  proposed  the  addi¬ 
tion  of  a  moving  ‘cap’  to  the  Drucker-Prager  model 
at  high-pressure  values.  The  cap  intersects  the 
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hydrostatic  loading  C h)  line  (to  account  for  the  fact  2.1.3  Flow  rule 
that  plastic  response  of  the  soil  can  be  attained  under 

pure  hydrostatic-pressure  loading  conditions).  The  ^he  third  component  of  the  material  model  is  the 

position  of  the  cap  is  taken  to  be  controlled  by  the  mag-  ^ow  ru*e’  w^'c^  sPcc'fies  ^1C  ev0'u  l'on  ^ie  inelastic 

nitude  of  the  plastic  volumetric  strain,  ep'  .  In  addition,  strain  components  during  deformation.  For  an  elastic- 

the  effect  of  negative  pressure  on  the  failure  of  soil  is  Plastic  material-model  formulation,  a  ‘normality  flow 

also  accounted  for  by  introducing  a  ‘tension  cut-off  rule js  generally  adopted,  which  can  be  defined  (in  a 

level  (the  level  of  negative  pressure  at  which  soil  fails  rate  orm'  as 

and  its  shear  strength  i.e.  drops  to  zero) .  The  pres-  .  g y 

sure  cut-off  is  denoted  as  T  in  Fig.  1(b).  Furthermore,  £pl  =  k  —  (9) 

to  account  for  the  fact  that,  as  pressure  is  increased, 

its  effect  on  shear  strength  decreases,  the  Drucker-  where  a  raised  dot  is  used  to  define  the  time  derivative, 

Prager  linear  (shear)  yield  function,  equation  (3),  is  superscript  pi  stands  for  plastic,  X  is  a  plastic  multi- 

modified  as  plier,  while  e  and  a  are  the  strain  and  Cauchy  stress 

(6  x  1)  column  vectors.  Equation  (9)  simply  states 
/shear C/2 >  h)  =  \fh  —  [a  —  y  exp(/i/j )  -  OR]  (4)  that  the  plastic  deformation  takes  place  in  a  direction 

which  is  normal  to  the  surface  defined  by  the  yield 
where  y  and  /3  are  additional  material  parameters.  function/ (i.e.  to  the  yield  surface). 

The  cap  portion  of  the  yield  function  is  assumed  to  Equation  (9)  can  be  rewritten  in  form  of  the  incre- 
be  elliptical  as  ments  of  the  plastic  strain  components,  Aepl ,  by  mul- 

1  tiplying  both  sides  of  equation  (9)  by  At  (a  small  time 

/cap  (jz,  h,  ep’^  =  Vh  -  —  y/(X  -  L)2  -  (/  -  L)2  (5)  increment).  For  a  rate-independent  (i.e.  inviscid)  plas¬ 

tic  material,  AX  =  X  At  is  constant  and  is  computed  (as 


where  R  is  a  material  parameter  (a  ratio  of  the  hori¬ 
zontal  and  the  vertical  axes  of  the  ellipse),  while 
X(ePgj)  and  are  the  plastic-volumetric-strain- 

dependent  intersection  of  the  cap  yield  function  and 
the  7i  axis  and  the  plastic  volumetric-strain  dependent 
initial  value  of  /  for  the  cap  yield-function. 


will  be  shown  in  next  section)  using  the  condition  that 
the  (new)  stress  state  (i.e.  the  stress  state  at  the  end 
of  an  elastic-plastic  loading  step  of  the  duration  At) 
satisfies  the  f(a)  —  0  condition.  In  the  case  of  visco¬ 
plastic  materials,  on  the  other  hand,  AX  is  not  constant 
and  is  usually  defined  as  [26] 


2. 1 .2  Constitutive  relation 


AX  —  rj(p(f)At 


(10) 


To  completely  define  the  position  of  the  cap  yield 
function  (which  changes  as  a  result  of  plastic  defor¬ 
mation),  a  plastic-volumetric-strain-dependent  hard¬ 
ening  parameter,  k  (the  value  of  /,  corresponding  to 
the  centre  of  the  cap  yield-function  ellipse)  is  used. 
This  parameter  is  defined  using  the  following  function 

x(4,)-k 

- J- - =a-y  exp(/i/i)  -  Of  (6) 

K 

while  the  relationship  between  k  and  L  is  given  as 

L  =  k,  if  k>  0 
L  =  0,  if  k  <  0 


where  r]  is  material  constant  commonly  referred  to  a 
fluidity  parameter  and  /(/)  is  defined  as 

**>-(£)"■  if^°  in, 

0(f)  =  0,  if/  ^  0 

where  f0  and  N  are  normalizing  constants  with  the 
units  of  /  and  an  exponent,  respectively.  Since,  in 
accordance  with  equation  (11),  AX  scales  with  At,  and 
there  is  no  requirement  that  the  new  stress  state  satis¬ 
fies  the  / (a)  —  0  condition,  the  inelastic  response  of 
the  corresponding  material  becomes  rate  (i.e.  time) 
dependent  (i.e.  the  material  become  visco-plastic). 


and  the  plastic  volumetric- strain  dependence  of  X  is 
furthermore  defined  as 


withD,  W,  and  Xn  =  X(ep*j  =  0)  being  material  param¬ 
eters  for  the  cap  yield-function.  It  should  be  noted  that 
the  portions  of  the  yield  surface  associated  with  shear 
and  tensile  cut-off  are  assumed  to  be  ‘ideal  plastic’,  i.e. 
not  to  evolve  during  plastic  deformation. 


2.2  Material-state  integration  procedure 

In  this  section,  a  brief  description  is  provided  of 
the  standard  procedure  that  is  used  to  evaluate  the 
response  of  an  elastic-plastic  or  an  elastic-visco- 
plastic  material  subjected  to  a  general  state  of  loading. 
Since  the  sand  model  developed  in  this  work  is  of 
an  elastic-visco-plastic  nature,  the  procedure  pro¬ 
vided  below  is  given  for  such  a  material  while,  where 
needed,  differences  with  the  elastic-plastic  behaviour 
are  emphasized. 
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First,  the  total  strain  increment  Ae  can  be  decom¬ 
posed  into  its  elastic  Aeel  and  visco-plastic  Ae'1’  as 

e  =  eel  4-  evp  (12) 

Next,  the  stress  at  the  end  of  an  elastic-visco-plastic 
deformation  step  (time  =  t  +  At)  is  defined  in  terms 
of  its  value  at  the  beginning  of  the  same  loading  step 
(time  =  t)  as 


the  associated  system  of  six  equations  for  an  elastic- 
visco-plastic  material,  equation  (13)  is  first  rewritten 
as 

C_1or(+At  +  A  tef+At  =  C~lat  +  Ae  (17) 

so  that  the  unknowns  are  located  on  the  left-hand  side, 
while  the  known  quantities  on  the  right-hand  side  of 
equation  (17).  Using  a  general  notation,  equation  (17) 
can  then  be  written  as 


ort+At  =<ft  +  C(Ae  -  Aevp) 


(13) 


where  C  is  the  6x6  elastic  stiffness  matrix  and 

Aevp  =  [(1  -  x)e?  +  X«^At]  Af  (14) 

where  0  <  x  <  1  is  an  adjustable  parameter  and  /  = 
0.0  corresponds  to  the  fully  explicit  formulation,  while 
X  =  1.0  corresponds  to  the  fully  implicit  formulation. 
The  value  of  x  controls  the  stability  of  the  numeri¬ 
cal  procedure  and  when  x  <  0.5  the  procedure  is 
conditionally  stable,  while  when  /  >  0.5  it  is  uncon¬ 
ditionally  stable  (i.e.  it  is  stable  for  any  value  of  the 
time  increment  At). 

In  the  fully- implicit  procedure  (x  =  1.0),  which  is 
considered  below,  the  visco-plastic  strain  increment  is 
defined  as 


Ae 


vp 

t+At 


AXt+At 


V 

3  a 


—  tl4>(ft+At)At 

t+At 


V 

da 


t+At 


(15) 


When  equation  (15)  is  plugged  into  equation  (13), 
a  system  of  six  non-linear  algebraic  equations  with 
six  unknown  components  of  <r(+A,  is  attained.  While 
one  of  the  standard  routines  can  be  invoked  to  solve 
this  equation  system,  a  computationally  more  efficient 
procedure  is  presented  below  [26] . 

When  a  material  is  elastic-plastic,  the  equation  (4) 
becomes 


Ae 


pi 

t+At 


epl 

* t+At 


At  — 


df 

AXt+At  - — 


da 


t+At 


(16) 


and  when  equation  (16)  is  plugged  into  equation  (13), 
a  system  of  six  equations  with  seven  unknowns  [at+At, 
AXt+At)  is  attained.  In  this  case,  the  seventh  equation 
is  obtained  by  setting/(<rf+At)  =  0. 


2.3  Numerical  solution  procedure 

In  this  section,  a  brief  description  is  provided  of  an 
efficient  numerical  procedure  that  is  used  to  inte¬ 
grate  the  response  of  an  elastic-visco-plastic  material 
subjected  to  a  general  state  of  loading  [26].  To  solve 


Pt+Ati^t+At)  —  <7 1 


(18) 


Since,  in  accordance  with  equation  (13),  evp  is  a  func¬ 
tion  of  stresses,  Pt+At  in  equation  (18)  is  implicitly  a 
function  of  at+At  alone. 

To  solve  equation  (18),  its  left-hand  side  is  next 
approximated  using  a  truncated  Taylor  series  as 

Pt+At  +  Pt+At  dat  +At  =  (19) 


where 


P 


f,i  _ 

t+At 


°rt+At 

d^t+At 


(20) 


and  i  (=  0, 1,2, . . .)  represents  the  iteration  number 
and  P';l+At  is  also  a  function  of  a\+At  alone. 

Thus,  when  an  initial  guess  is  made  for  a\+AV  then 
equation  (20)  can  be  solved  for  d<rj+At  and  an  updated 
value  for  stress,  a\  pAf,  obtained  as 


oftit  =  ®\+At  +  dff‘t+At 


Next,  the  total  visco -plastic  strain  el'^At 


e 


vp,/+l 

t+At 


—  et+At  C  ® 


i+ 1 
t+At 


(21) 

is  computed  as 
(22) 


and  since  the  yield  function  depends  on  both  the 
current  values  of  stress  state  and  the  total  visco¬ 
plastic  strain  (i.e.  the  stress  state  controls  the  values 
of  stress  invariants  while  the  total  visco -plastic  strain 
controls  the  hardening- dependent  cap  position),  it  is 
updated  as 

ft%=f(< Iv^a')  (23) 

Similarly,  equation  (15)  can  be  used  to  obtain  gf p 1 
and  since  both  oj ;  ]AI  and  e^lAtl  are  now  known,  P‘  [ 
and  P';[  A/  can  be  updated  and  the  aforementioned  pro¬ 
cedure  continued  until  a  convergence  |d<rj+Af  |  <  tol  is 
attained. 

It  should  be  noted  that,  as  correctly  pointed  out 
by  one  of  the  reviewers  of  the  present  article,  the 
solution  method  employed  in  the  present  work  is 
a  Newton-Raphson  stress-based  iteration  scheme. 
There  are,  however,  other  solution  procedures  such  as 
the  so-called  ‘cutting  plane’  purely  explicit  algorithm 
that  is  based  on  the  plastic  coefficient  lambda  and 
other  implicit  schemes  based  on  either  stress  or 
hardening  parameters  [27], 
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2.4  Derivation  and  parameterization  of  the 
CU-ARL  visco-plastic  sand  model 

In  this  section,  a  procedure  is  presented  for  the  exten¬ 
sion  of  the  elastic-visco-plastic  material  model  for 
sand  recently  proposed  by  Tong  and  Tuan  [26]  in 
order  to  include  the  effect  of  saturation  level  on  the 
sand  behaviour.  A  summary  of  the  proposed  mathe¬ 
matical  expressions  that  describe  the  effect  of  degree 
of  saturation  on  various  parameters  of  the  elastic- 
visco-plastic  material  model  for  sand  is  provided  in 
Table  1.  These  expressions  are  parameterized  using 
open-literature  experimental  data  and  a  property- 
correlation  analysis  and  the  parameterization  results 
are  also  given  in  Table  1 .  Derivation  and  parameteri¬ 
zation  are  carried  out  separately  for  the  yield  function, 
the  constitutive  equation,  and  the  visco -plastic  flow 
rule. 

2.4. 1  Isotropic  linear  elastic  response 

Since  the  shear  modulus  of  sand  is  mainly  con¬ 
trolled  by  the  shear  stiffness  of  its  skeleton  formed 
by  the  interlocked  particles,  the  degree  of  satura¬ 
tion  is  not  expected  to  have  a  first-order  effect 
on  the  magnitude  of  the  shear  modulus.  On  the 
other  hand,  since  the  pore-pressure  in  saturated 
sand  provides  additional  hydrostatic  pressure  load¬ 
carrying  capacity  to  the  sand,  the  bulk  modulus  is 
expected  to  increase  with  an  increase  in  the  degree 
of  saturation  at  both  the  low  and  high  deformation- 
rates. 


2.4.2  Yield  criterion 

As  seen  in  Fig.  1(b),  the  yield  surface  consists  of 
three  separate  parts:  the  shear  yield  surface,  the  cap 
yield  surface,  and  the  tensile  cut-off  limit.  These  are 
considered,  one  at  a  time,  in  this  section. 

Shear-yield  surface:  according  to  equation  (4),  this 
portion  of  the  yield  function  contains  four  parame¬ 
ters:  a,  p,  y,  andd.  a-y  is  the  zero-pressure  shear-yield 
strength  of  the  material,  y  and  fi  control  the  rate  of  loss 
of  the  pressure  dependence  of  the  shear-yield  strength 
while  6  quantifies  pressure  dependence  of  the  shear- 
yield  strength.  For  dry  sand,  these  parameters  were 
determined  by  Tong  and  Tuan  [26]  using  available 
public-domain  experimental  results,  e.g.  references 
[28]  and  [29].  These  parameters  for  dry  sand  were 
first  verified  in  the  present  work.  Next,  the  effect  of 
water /moisture  on  the  value  of  these  parameters  is 
considered. 

Since  water  present  in  sand  increases  adhesion 
between  the  sand  particles,  the  zero-pressure  shear 
strength  {—a-y)  is  expected  to  increase  with  an 
increase  in  the  degree  of  saturation.  However,  since 
water  also  acts  as  an  inter-particle  lubricant,  it  reduces 
the  inter-particle  friction  coefficient  and,  in  turn, 
the  pressure  dependence  of  the  shear  strength  (=  0). 
Lastly,  p  (and  y)  are  adjusted  to  comply  with  the  con¬ 
dition  that  there  is  a  smooth  transition  between  the 
shear  yield  and  the  cap  yield  functions. 

Cap-yield  function:  the  cap-yield  function, 
equation  (5),  is  essentially  controlled  by  two  param¬ 
eters:  X  (or  L )  and  R.  X  represents  the  (pure)  pressure 


Table  1  Parameterization  of  the  CU-ARL  visco-plastic  sand  model  for  sand  with  an 
initial  porosity  of  0.36 


Parameter 

symbol 

Unit 

Equations  where 
first  used 

Value/expression 
(dry  sand) 

value/expression 
(saturated  sand) 

Yield  function 

a 

MPa 

(4) 

0.0642 

0.0642(1  -  ftat)  +  0.077/W 

Y 

MPa 

(4) 

0.00  589 

0.00  589(1  -  ftat)  +  0.0062fta, 

P 

MPa-1 

(4) 

0.34  283 

0.34283 

e 

- 

(4) 

0.18  257 

0.18257(1 -fta,)  +  0.1187ftat 

T 

MPa 

0.0069 

0.0069(1  -  ftat)  +  0.0084ftat 

Constitutive  relation 

W 

(8) 

0.2142 

0.2142 

D 

MPa-1 

(8) 

0.00  952 

0.00  952 

R 

- 

(6) 

5.0 

5.0 

X0 

MPa 

(8) 

0.01 

0.01(1  -  ftat)  +  0.0121fta, 

Flow  rule 

^low-rate 

|TS_1 

(10) 

0.0002 

0.0002(1  -ftat)  +  0.00  016fta, 

*7  high -rate 

(TS-1 

(10) 

- 

0.002(1  —  ftat)  +  0.0016ftat 

N 

- 

(11) 

1.0 

1.0 

fo 

MPa 

(11) 

100 000 

100  000 

Isotropic  linear  elastic  response 

G  MPa  (13) 

63.85 

63.85(1  -  ft at)  +  61.28ftat 

^low-rate 

MPa 

(13) 

106.4 

106.4(1  -  ftaf)  +  127.68ftat 

^high-rate 

MPa 

(13) 

3000.0 

3000.0(1  -  ftat)  +  3900.00ftat 

ftat  =  degree  of  saturation. 
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under  which  plastic  compaction  takes  place.  Since 
plastic  compaction  relies  on  the  motion  of  sand 
particles,  and  thus  is  affected  by  inter-particle  adhe¬ 
sion/friction,  the  initial  value  of  X,  Xn,  is  expected  to 
increase  with  the  corresponding  zero-pressure  value 
of  the  shear-yield  strength.  As  far  as  the  cap-yield 
ellipse  horizontal-to-vertical  axis  ratio  (eccentricity), 
R,  is  concerned,  it  is  assumed  that  it  does  not  change 
with  the  degree  of  saturation  of  the  sand. 

Tension  cup-off:  the  tensile  fracture  of  sand  involves 
inter-particle  bond  breaking  and,  hence,  the  tensile 
cut-off  limit,  T,  is  expected  to  be  governed  by  inter¬ 
particle  adhesion  and  scale  with  the  zero-pressure 
shear-yield  strength. 

2.4.3  Constitutive  relation 

The  basic  assumption  used  by  Tong  and  Tuan  [26] 
that  only  the  cap  position  of  the  yield  surface  evolves 
during  plastic  deformation  and  that  only  the  vol¬ 
umetric  part  of  the  plastic  deformation  causes  the 
evolution  of  the  yield  surface,  is  retained.  According 
to  equation  (8),  the  constitutive  relation  for  sand  con¬ 
tains  two  parameters  D  and  W  (X0  was  defined  by 
the  cap  yield  function) .  Limited  experimental  results 
available  for  parameterization  of  the  CU-ARL  visco¬ 
plastic  sand  model  did  not  allow  reliable  evaluation 
of  the  effect  of  saturation  on  D  and  W.  Conse¬ 
quently,  in  the  present  rendition  of  the  model,  D  and 
W  are  assumed  to  be  independent  of  the  degree  of 
saturation. 

2.4.4  Flow  rule 

According  to  equations  (10)  and  (11),  the  visco-plastic 
flow  rule  function  contains  three  parameters:  i],  f0, 
and  N.  These  parameters  essentially  control  the  rate 
at  which  visco-plastic  deformation  takes  place  as  a 
function  of  the  driving  force  (/  >  0),  that  is:  r]/f™ 
acts  as  a  mobility  term  while  N  is  a  viscous-exponent 
factor.  Since  water  acts  as  an  inter-particle  lubricant 
and,  at  the  same  time,  as  an  inter-particle  cohesion 
agent,  the  effect  of  water  in  sand  on  these  param¬ 
eters  is  not  easy  to  predict.  Also,  the  effect  can  be 
different  at  different  deformation  rates,  since  at  high 
deformation  rates  water  is  trapped  and  less  volumet¬ 
ric  plastic  deformation  takes  place.  On  the  other  hand, 
at  lower  deformation  rates,  water  is  typically  squeezed 
out  and  more  plastic  compaction  is  observed.  Limited 
experimental  data  available  for  model  parameteriza¬ 
tion  allowed  only  a  rough  assessment  of  the  effect 
of  the  degree  of  saturation  on  the  flow  rule  relation. 
Within  this  assessment,  f0  and  N  were  found  not  to 
be  significantly  affected  by  the  degree  of  saturation. 
Also,  i]  was  found  not  to  be  significantly  affected 
by  the  degree  of  saturation  but  only  in  the  high- 
deformation-rate  regime.  In  the  low-deformation-rate 
regime,  ?;  was  found  to  increase  significantly  with 


the  increase  in  degree  of  saturation  suggesting  that 
saturated  sand  tends  to  behave  more  like  a  plastic 
material. 

2.4.5  Model  parameterization 

The  procedure  presented  above  was  applied  to  the 
case  of  sand  with  an  initial  porosity  of  0.36.  The  result¬ 
ing  parameterization  of  the  CU-ARL  visco-plastic  sand 
model  is  given  in  Table  1. 

2.5  Implementation  of  the  CU-ARL  visco-plastic 
sand  model 

The  CU-ARL  visco-plastic  sand  model  described  in 
the  previous  section  is  next  implemented  in  the 
material  user  subroutine,  VUMAT,  of  the  commer¬ 
cial  finite -element  program  ABAQUS/ Explicit  [30], 
This  subroutine  is  compiled  and  linked  with  the 
finite-element  solver  and  enables  ABAQUS /Explicit  to 
obtain  the  needed  information  regarding  the  state  of 
the  material  and  the  material  mechanical  response 
during  each  time  step,  for  each  integration  point 
of  each  element.  In  the  present  work  first-order 
eight-node  general-purpose  reduced-integration  solid 
elements  (ABAQUS /Explicit  designation  C3D8R)  are 
used. 

The  essential  features  of  the  coupling  between 
the  ABAQUS /Explicit  finite- element  solver  and  the 
VUMAT  Material  User  Subroutine  at  each  time  incre¬ 
ment  at  each  integration  point  of  each  element  can  be 
summarized  as  follows. 

1.  The  corresponding  previous  time-increment 
stresses  and  material  state  variables  as  well  as  the 
current  time-step  deformation  gradient  are  pro¬ 
vided  by  the  ABAQUS /Explicit  finite- element  solver 
to  the  material  subroutine.  In  the  present  work 
six  total  strain  components,  the  volumetric  plas¬ 
tic  strain,  six  visco -plastic  strain-rate  components, 
and  an  element-deletion  status  flag  were  used  as 
the  state  variables. 

2.  Using  the  information  provided  in  step  1,  and  the 
CU-ARL  visco-plastic  sand  model  presented  in  the 
previous  section,  the  material  stress  state  as  well 
as  values  of  the  material  state  variables  at  the 
end  of  the  time  increment  are  determined  within 
the  VUMAT  and  returned  to  the  ABAQUS /Explicit 
finite-element  solver.  In  addition,  the  changes  in 
the  total  internal  and  the  inelastic  energies  (where 
appropriate)  are  computed  and  returned  to  the 
solver. 

It  should  be  noted,  that  the  VUMAT  file  developed 
in  the  present  work  could  be  used  with  minor  mod¬ 
ifications  and  new  parameterization  for  a  variety  of 
soils  with  different  chemical  make-up  and  granular 
microstructure. 


JMDA237  ©  IMechE  2009 


Proc.  IMechE  Vol.  223  Part  L:  J.  Materials:  Design  and  Applications 


70 


M  Grujicic,  T  He,  B  Pandurangan,  W  C  Bell,  B  A  Cheeseman,  W  N  Roy,  and  R  R  Skaggs 


3  VALIDATION  OF  THE  CU-ARL  VISCO¬ 
PLASTIC  SAND  MODEL  IN  THE 
HIGH-DEFORMATION-RATE  REGIME 

In  this  section,  the  CU-ARL  visco-plastic  sand  model 
is  tested  by  carrying  out  a  number  of  computational 
transient  non-linear  dynamics  analyses  under  high 
deformation-rate  conditions  and  by  comparing  the 
computed  results  with  their  experimental  counter¬ 
parts.  All  the  calculations  carried  out  in  this  section 
were  done  using  ABAQUS/ Explicit  [30],  a  general- 
purpose  non-linear  dynamics  modelling  and  simula¬ 
tion  software.  In  this  section,  a  brief  overview  is  given 
of  the  basic  features  of  ABAQUS /Explicit,  emphasizing 
the  aspects  of  this  computer  program  that  pertain  to 
the  problem  at  hand. 

A  transient  non-linear  dynamics  problem  is  anal¬ 
ysed  within  ABAQUS /Explicit  by  solving  simultane¬ 
ously  the  governing  partial  differential  equations  for 
the  conservation  of  momentum,  mass,  and  energy 
along  with  the  materials  constitutive  equations  and 
the  equations  defining  the  initial  and  the  bound¬ 
ary  conditions.  The  equations  mentioned  above  are 
solved  numerically  using  a  second-order  accurate 
explicit  scheme  and  one  of  the  two  basic  mathe¬ 
matical  approaches,  the  Lagrange  approach  and  the 
Euler  approach.  The  key  difference  between  the  two 
approaches  is  that  within  the  Lagrange  approach 
the  numerical  grid  is  attached  to  and  moves  along 
with  the  material  during  calculation,  while  within 
the  Euler  approach,  the  numerical  grid  is  fixed  in 
space  and  the  material  moves  through  it.  Within 
ABAQUS /Explicit,  the  Lagrange  approach  is  used. 
In  our  recent  work  [9],  a  brief  discussion  was 
given  of  how  the  governing  differential  equations 
and  the  materials  constitutive  models  define  a  self- 
consistent  system  of  equations  for  the  dependent 
variables  (nodal  displacements,  nodal  velocities,  ele¬ 
ment  material  densities,  and  element  internal  energy 
densities) . 

Throughout  this  article,  the  terms  ‘depth  of  burial’ 
(DOB)  and  ‘stand-off  distance’  (SOD)  are  used  to 
denote  distances  between  the  mine  top-face  and 
the  sand  top-face  and  between  the  sand  top- 
face  and  the  bottom  face  of  the  target  structure, 
respectively. 


3. 1  Total  momentum  transferred  to  the 
target  structure 

To  assess  the  ability  of  the  CU-ARL  visco-plastic  sand 
model  to  account  for  the  total  momentum  trans¬ 
ferred  to  the  target  structure  following  detonation  of 
a  ground-laid  or  shallow-buried  mine  at  different  sat¬ 
uration  levels  of  the  sand,  the  computational  results 
are  compared  with  their  experimental  counterparts 
obtained  in  references  [8]  and  [31]. 


3.1.1  Dry  sand 

To  assess  the  ability  of  the  CU-ARL  visco-plastic 
sand  model  to  account  for  the  total  momentum 
transferred  to  the  target  structure  at  low  satura¬ 
tion  levels  of  the  sand,  a  non-linear  dynamics-based 
computational  analysis  of  the  interaction  of  detona¬ 
tion  products,  mine  fragments,  and  sand  ejecta  with 
an  instrumented  horizontal  mine-impulse  pendulum 
used  in  reference  [8]  is  carried  out  and  the  computed 
results  compared  with  their  experimental  counter¬ 
parts.  Since  an  in-depth  description  of  the  experimen¬ 
tal  details  related  to  the  construction  and  utilization  of 
the  instrumented  horizontal  mine-impulse  pendulum 
can  be  found  in  our  recent  work  [32] ,  they  will  not  be 
presented  here. 

Next,  a  brief  description  is  given  of  the  compu¬ 
tational  model  used  to  simulate  the  interaction  of 
the  detonation-products/soil  ejecta  resulting  from  the 
explosion  of  a  shallow-buried  or  ground-laid  mine  and 
the  instrumented  horizontal  mine-impulse  pendu¬ 
lum.  The  computational  modelling  of  this  interaction 
involves  two  distinct  steps:  (a)  geometrical  modelling 
of  the  instrumented  horizontal  mine-impulse  pendu¬ 
lum  set-up  and  (b)  a  non-linear  dynamics  analysis 
of  the  momentum  transfer  from  the  detonation- 
products,  the  mine  casing,  and  soil  ejecta  to  the 
pendulum. 

The  three  (pendulum,  sand,  and  mine)  computa¬ 
tional  domains  used  in  the  present  study  are  shown 
in  Fig.  2.  The  geometrical  models  for  the  various 
components  of  the  pendulum  were  constructed  using 
eight-node  reduced-integration  solid  elements  with 
sizes  varying  between  50  x  50  x  20  mm  and  100  x 
50  x  50  mm.  An  advantage  was  taken  of  the  planar 
symmetry  of  the  model.  In  other  words,  a  vertical  plane 
of  symmetry  was  placed  along  the  length  of  the  pen¬ 
dulum  that  enabled  only  a  half  of  the  pendulum,  the 
sand,  and  the  mine  to  be  modelled.  In  accordance 
with  the  instrumented  horizontal  mine-impulse  pen¬ 
dulum  used  in  reference  [8] ,  different  sections  of  the 
pendulum  were  constructed  using  AISI  1006  steel 
and  rolled  homogenized  armour  (RHA)  plate  material. 
Welded  joints  of  the  different  sections  of  the  pendu¬ 
lum  were  simulated  by  joining  the  components  in 
question. 

The  size  and  circular-disc  shape  of  the  mine  com¬ 
putational  domain  were  selected  to  match  the  ones 
for  the  C4  mine  used  in  reference  [8].  The  mine 
computational  domain  was  meshed  using  eight-node 
reduced-integration  solid  elements  with  a  typical  size 
of  5  x  5  x  5  mm  and  filled  with  a  C4  HE  material. 

The  sand  computational  domain  was  modelled  as 
a  circular  cylinder  with  a  radius  of  750  mm  and  a 
height  of  400  mm.  The  domain  was  divided  into  three 
concentric  subdomains.  All  three  subdomains  were 
meshed  using  eight-node  reduced-integration  solid 
elements  with  a  typical  mesh  size  of  5  x  5  x  10  mm  in 
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Fig.  2  Computational  domains  used  in  the  analysis  of  total  momentum  transferred 
to  the  horizontal  mine-impulse  pendulum  in  the  case  of  dry  sand 


the  innermost  subdomain  and  a  maximum  mesh  size 
of  50  x  50  x  iO  mm  in  the  outer-most  sub-domain. 
The  lateral  and  the  bottom  faces  of  the  sand  domain 
were  lastly  surrounded  with  eight-node  CIN3D8  infi¬ 
nite  elements  in  order  to  model  far-held  sand  regions 
and  avoid  un-physical  stress-wave  reflection  at  the 
sand-domain  lateral  and  bottom  surfaces.  The  sand 
domains  containing  C3D8R  elements  were  filled  with 
CU-ARL  visco-plastic  sand  material  while  the  infinite 
elements  were  filled  with  an  ‘elastic’  sand  material 
with  a  Young’s  modulus  and  a  Poisson’s  ratio  match¬ 
ing  those  of  the  CU-ARL  visco-plastic  sand.  Also,  to 
account  for  non-linear  behaviour  of  the  bulk  mod¬ 
ulus  at  high  deformation  rates  and  large  volumetric 
strains,  and  in  accordance  with  the  procedure  used  in 
the  Laine  and  Sandvik  model  [5] ,  a  higher  value  (3  GPa, 
Table  1)  of  the  bulk  modulus  was  used  in  this  portion 
of  the  work. 

The  mine-sand  and  sand-pendulum  interactions 
were  modelled  using  the  ‘hard  contact  pair’  type  of 
contact  algorithm.  Within  this  algorithm,  contact  pres¬ 
sures  between  two  bodies  are  not  transmitted  unless 
the  nodes  on  the  ‘slave  surface’  contact  the  ‘master 
surface’.  No  penetration /over  closure  is  allowed  and 
there  is  no  limit  to  the  magnitude  of  the  contact  pres¬ 
sure  that  could  be  transmitted  when  the  surfaces  are 
in  contact.  Transmission  of  shear  stresses  across  the 
contact  interfaces  is  defined  in  terms  of  a  static  and 
a  kinematic  friction  coefficient  and  an  upper-bound 
shear  stress  limit  (a  maximum  value  of  shear  stress 
that  can  be  transmitted  before  the  contacting  surfaces 
begin  to  slide) . 

A  standard  mesh- sensitivity  analysis  was  carried  out 
(the  results  not  shown  for  brevity)  in  order  to  ensure 
that  the  results  obtained  are  insensitive  to  the  size  of 
the  elements  used.  Similar  mesh-sensitivity  analyses 
were  carried  out  for  the  remaining  studies  presented 
in  this  article. 


At  the  beginning  of  the  simulation,  the  pendulum 
was  assumed  to  be  at  rest  (with  the  gravitational  force 
acting  downwards),  while  the  sand  and  mine  domains 
were  filled  with  stationary  materials  (sand  and  C4, 
respectively).  Mine  detonation  was  initiated  first  along 
the  bottom  face  of  the  mine.  The  motion  of  the  pendu¬ 
lum  was  constrained  to  within  a  vertical  plane  and  a 
fixed  single-point  constraint  was  applied  to  its  pivot 
point.  The  fixed  boundary  conditions  were  applied 
along  the  circumferential  and  the  bottom  faces  of  the 
sand  domain. 

fn  accordance  with  the  experimental  procedure 
employed  in  reference  [8],  the  resultant  maximum 
angular  displacement  of  the  pendulum  arm  is  mea¬ 
sured  and  used  to  calculate  the  detonation-induced 
impulse  on  the  pendulum.  More  details  pertaining  to 
the  procedure  used  to  assess  the  total  mine  impulse 
transferred  to  the  pendulum  can  be  found  in  our 
previous  work  [32] . 

An  example  of  the  results  pertaining  to  the  tem¬ 
poral  evolution  of  the  materials  distribution  in  the 
pendulum,  the  sand  domain  and  the  mine  domain 
obtained  in  the  present  work,  is  displayed  in  Figs  3(a) 
and  (b).  Clearly,  upon  detonation  of  the  mine,  a  sand- 
overburden  bubble  is  formed  that  expands  upward 
and  laterally.  When  the  bubble  makes  contact  with 
the  RHA  and  the  impact  plates  of  the  pendulum, 
the  momentum  transfer  from  sand  to  the  pendulum 
begins  and  the  pendulum  swings  upward. 

The  effect  of  the  DOB  of  the  mine  on  the  total 
momentum  transferred  to  the  pendulum  in  the  case 
of  dry  sand  obtained  in  reference  [32]  is  displayed  in 
Fig.  4.  Also  displayed  in  Fig.  4  are  the  computational 
results  obtained  in  the  present  work  using  both  the 
original  sand  compaction  model  of  Laine  and  Sandvik 
[5]  and  the  present  CU-ARL  visco-plastic  sand  model. 
The  results  displayed  in  Fig.  4  can  be  summarized  as 
follows. 
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Fig.  3  An  example  of  the  temporal  evolution  of  mate¬ 
rials’  distribution  during  the  computational 
analysis  of  total  momentum  transferred  to  the 
horizontal  mine-impulse  pendulum  in  the  case 
of  dry  sand.  Post  detonation  times:  (a)  230  ps  and 
(b)  375  ps.  DOB  =  50  mm 


Fig.  4  The  effect  of  mine  DOB  on  the  total  impulse  trans¬ 
ferred  to  the  horizontal  mine-impulse  pendulum 
in  the  case  of  dry  sand 

1.  The  computation/experiment  agreement  is  some¬ 
what  improved  when  the  original  porous-material 
model  is  replaced  with  the  CU-ARL  visco-plastic 
sand  model  for  all  landmine  detonation  cases 
analysed. 


2.  Both  sets  of  computational  results,  however,  under¬ 
predict  the  total  momentum  transferred  to  the 
pendulum  at  both  DOBs.  This  finding  is  not  that 
unexpected  since  it  is  generally  observed  (e.g.  ref¬ 
erence  [32])  that  about  20-30  per  cent  of  the  total 
momentum  is  transferred  to  the  target  structure 
by  blast  waves  generated  at  the  sand-air  interface. 
Since,  in  the  present  work,  air  surrounding  the  pen¬ 
dulum  was  not  modelled,  one  could  expect  that  the 
computed  results  pertaining  to  the  total  momen¬ 
tum  transfer  are  lower  than  their  experimental 
counterparts  by  about  20-30  per  cent.  This  is  in  fact 
the  level  of  computation/ experiment  disagreement 
observed  in  Fig.  4. 

3.  The  general  increase  in  the  total  momentum  trans¬ 
ferred  to  the  target  structure  with  an  increase  in  the 
DOB  is  fairly  well  predicted  by  both  sand  models. 

4.  The  effect  of  the  DOB  on  the  momentum  transfer  is 
directly  related  to  the  extent  of  sand  ejection  and  to 
the  energy  absorbed  by  the  sand  surrounding  the 
mine.  Since  rate  of  increase  of  the  total  momen¬ 
tum  transferred  to  the  pendulum  decreases  with 
an  increase  in  the  DOB,  it  appears  that  there  is  an 
optimum  DOB  that  maximizes  the  lethal  effect  of 
detonation  of  a  shallow-buried  mine.  This  can  be 
rationalized  by  the  fact  that  as  the  DOB  is  increased 
the  effects  of  detonation  become  more  confined 
within  the  soil  (the  ‘camouflet  effect’) .  Determina¬ 
tion  of  the  optimal  DOB  was  beyond  the  scope  of 
the  present  work. 


3.1.2  Saturated  sand 

To  assess  the  ability  of  the  CU-ARL  visco -plastic  sand 
model  to  account  for  the  total  momentum  trans¬ 
ferred  to  the  target  structure  at  high  saturation  levels 
of  the  sand,  a  non-linear  dynamics-based  computa¬ 
tional  analysis  of  the  interaction  of  detonation  prod¬ 
ucts,  mine  fragments,  and  sand  ejecta  with  a  vertical 
impulse  measurement  fixture  (VIMF)  used  in  refer¬ 
ence  [31]  is  carried  out  and  the  computed  results 
compared  with  their  experimental  counterparts.  Since 
a  detailed  description  of  the  experimental  procedure 
related  to  the  construction  and  utilization  of  the  VIMF 
can  be  found  in  our  recent  work  [9] ,  they  will  not  be 
presented  here.  Rather,  it  could  be  stated  that  the  main 
component  of  the  VIMF  is  a  witness  plate  placed  at 
a  SOD  from  the  sand  right  above  the  mine  burial- 
site.  Upon  detonation  of  the  mine,  the  interactions 
between  the  sand-overburden  bubble  and  the  plate 
results  in  momentum  transfer  that  is  assessed  from 
the  maximum  upward  velocity/ displacement  of  the 
plate  and  the  plate’s  mass. 

The  basic  formulation  of  the  computational  prob¬ 
lem  dealing  with  the  interactions  between  the  det¬ 
onation  products,  shell  fragments,  and  soil  ejecta 
(all  resulting  from  the  explosion  of  a  shallow-buried 
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landmine)  and  the  VIMF  is  presented  next.  The  com¬ 
putational  modelling  of  this  interaction  involved  two 
distinct  steps:  (a)  geometrical  modelling  of  the  VIMF 
along  with  the  adjoining  mine  and  sand  regions  and 
(b)  the  associated  transient  non-linear  dynamics  ana¬ 
lysis  of  the  impulse  loading  (momentum  transfer) 
from  the  detonation  products,  shell  fragments,  and 
soil  ejecta  to  the  VIMF  structure. 

The  three  (VIMF  witness  plate,  sand,  and  mine) 
computational  domains  used  in  the  present  study  are 
shown  in  Fig.  5.  Due  to  inherent  symmetry  of  the 
model,  the  problem  is  treated  as  being  axisymmet- 
ric.  The  geometrical  model  for  the  VIMF  witness  plate, 
sand,  and  mine  were  all  constructed  using  four-node 
axisymmetric  CAX4  solid  elements  with  sizes  varying 
between  2x2  mm  and  30  x  30  mm.  The  lateral  and 
the  bottom  faces  of  the  sand  domain  were  lastly  sur¬ 
rounded  with  four-node  CINAX4  infinite  elements  in 
order  to  model  far-held  sand  regions  and  avoid  un¬ 
physical  stress-wave  reflection  at  the  sand-domain 
lateral  and  bottom  surfaces.  The  VIMF  witness  plate 
subdomain  was  filled  with  AISI  4340  steel,  the  sand 
subdomain  was  filled  with  the  present  CU-ARL  visco¬ 
plastic  sand  material  or  its  elastic  equivalent,  and  the 
mine  subdomain  was  filled  with  TNT  HE  material.  The 
VIMF  and  the  mine  material  models  and  material- 
model  parameters  could  be  found  in  our  previous 
work  [9]. 

The  mine-sand  and  sand-VIMF  witness  plate  inter¬ 
actions  were  modelled  using  the  aforementioned  hard 
contact  algorithm.  The  total  mine  impulse  transferred 
to  the  VIMF  witness  plate  was  calculated  using  the 
maximum  average  vertical  displacement  of  the  plate. 
The  remainder  of  the  computational  procedure  fol¬ 
lowed  the  steps  outlined  previously  in  the  conjunc¬ 
tion  with  the  mine-impulse  pendulum  momentum 
transfer  analysis. 

An  example  of  the  computation/experiment  com¬ 
parison  carried  out  in  this  portion  of  the  work  is 
displayed  in  Fig.  6  in  which  the  effect  of  DOB  on 
the  total  momentum  transferred  to  the  VIMF  witness 
plate  is  displayed  for  the  case  of  a  4.54  kg,  circular-disc 


Fig.  5  Computational  domains  used  in  the  analysis  of 
total  momentum  transferred  to  the  VIMF  witness 
plate  in  the  case  of  saturated  sand 


Fig.  6  The  effect  of  mine  DOB  on  the  total  impulse  trans¬ 
ferred  to  the  VIMF  witness  plate  in  the  case  of 
saturated  sand 


TNT  mine  and  SOD  =  400  mm.  A  brief  analysis  of  the 
results  displayed  in  Fig.  6  reveals  the  following. 

1.  A  significant  improvement  in  the  computa¬ 
tion/experiment  agreement  is  obtained  when  the 
compaction  model  of  Laine  and  Sandvik  [5]  is 
replaced  by  the  CU-ARL  visco-plastic  sand  model. 
This  finding  was  expected  since  the  compaction 
model  does  not  include  the  effect  of  moisture  in 
sand.  This  is  found  to  result  in  excessive  energy 
absorption  by  the  dry  sand  due  to  its  compaction 
and  in  a  significantly  reduced  extent  of  sand  ejec¬ 
tion  (both  of  these  effects  result  in  reduced  levels  of 
total  momentum  transfer  to  the  target  structure). 

2.  The  total  momentum  transfer  values  predicted  by 
the  CU-ARL  visco-plastic  sand  model  are  again 
lower  by  ~ 20-30  per  cent  than  their  experimental 
counterparts.  As  stated  earlier,  this  discrepancy  is 
related  to  lack  of  consideration  of  an  air  domain 
connecting  the  sand  domain  and  the  pendulum, 
and  due  to  the  associated  neglect  of  the  air-borne 
blast  waves  generated  at  the  sand-air  interface  that 
impact  the  VIMF  witness  plate  before  sand  ejecta 
reach  the  same  plate. 

3.2  Spatial  and  temporal  evolution  of  the 
sand-overburden  bubble 

To  further  assess  the  validity  of  the  CU-ARL  visco¬ 
plastic  sand  model,  spatial  and  temporal  evolutions  of 
the  sand- overburden  bubble  (i.e.  the  elevated  portion 
of  sand  located  above  the  mine)  following  detona¬ 
tion  of  a  shallow  buried  mine  in  dry  and  saturated 
sand  were  analysed  computationally  and  the  com¬ 
putational  results  compared  with  their  experimen¬ 
tal  counterparts  obtained  in  reference  [33].  In  this 
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section,  a  brief  overview  of  the  experimental  setup  and 
the  procedure  used  in  reference  [33]  is  first  presented. 

The  experiments  carried  out  in  reference  [33]  can  be 
briefly  described  as  follows:  a  12.7  mm  wall  thickness 
cylindrical  barrel  with  an  outer  diameter  of  810  mm 
and  an  overall  height  of  710  mm  is  filled  with  sand  up 
to  its  top.  A 1 00  g  cylindrical- disc  shape  C4  high- energy 
explosive  (64  mm  in  diameter  and  20  mm  in  height)  is 
buried  into  the  sand  along  the  centre-line  of  the  bar¬ 
rel  with  its  faces  parallel  with  the  sand  surface.  The 
DOB’s  of  30  mm  and  80  mm  were  investigated.  The 
temporal  evolutions  of  the  diameter  and  the  height 
of  the  sand-overburden  bubble  were  determined  by 
analysing  the  data  obtained  using  high-speed  X-ray 
photography. 

The  computational  model  developed  in  this  portion 
of  the  work  included  three  subdomains:  (a)  a  steel- 
barrel  subdomain,  (b)  a  sand  subdomain,  and  (c)  a 
C4  mine  subdomain.  Again,  the  problem  was  con¬ 
sidered  as  being  axisymmetric.  The  steel  barrel  was 
modelled  using  four-node  axisymmetric  CAX4  solid 
elements  with  a  typical  size  of  6.35  x  5  mm  and  filled 
with  AISI  4340  steel.  No  boundary  conditions  were 
prescribed  along  the  outer  surfaces  of  the  steel  bar¬ 
rel.  The  remainder  of  the  computational  domain  was 
constructed  using  a  procedure  analogous  to  the  one 
previously  described  for  the  VIMF  setup.  Due  to  the 
similarity  of  the  present  computational  model  and  the 
one  used  for  the  VIMF  analysis,  no  further  descrip¬ 
tion  of  the  present  model  will  be  provided  and  no 
figures  dealing  with  the  finite -element  mesh  used  will 
be  displayed. 

An  example  of  the  typical  effect  of  the  presence  of 
moisture  on  the  morphology  of  the  sand-overburden 
bubble  is  displayed  in  Figs  7(a)  and  (b).  In  excellent 
agreement  with  the  general  experimental  observa¬ 
tions,  the  sand  bubble  is  shorter  and  wider  in  the  case 
of  dry  sand,  Fig.  7  (a) ,  and  taller  and  narrower,  Fig.  7  (b) , 
in  the  case  of  saturated  sand. 

3.2. 1  Dry  sand 

The  ability  of  the  CU-ARL  visco-plastic  sand  model 
to  account  for  the  main  observations  obtained  during 
experimental  investigation  of  landmine  detonation  in 
‘dry’  sand  (average  degree  of  saturation  /Jsat  <  0.15) 
[33]  is  discussed  first.  A  comparison  of  the  compu¬ 
tational  and  experimental  results  pertaining  to  the 
effect  of  SOD  on  the  temporal  evolution  of  the  dry 
sand-overburden  bubble  height  and  radius  is  shown 
respectively  in  Figs  8(a)  and  (b).  The  results  displayed 
in  Figs  8(a)  and  (b)  show  that  the  overall  agreement 
between  the  CU-ARL  visco-plastic  sand  model-based 
results  and  the  experimental  results  is  reasonably  good 
particularly  in  the  case  of  30  mm  DOB.  The  results 
obtained  using  the  compaction  model  showed  signif¬ 
icantly  worse  agreement  with  the  experiment  and,  to 
keep  clarity  of  the  figures,  are  not  displayed  in  Figs  8(a) 


Fig.  7  Typical  morphology  of  the  sand-overburden  bub¬ 
ble  in  the  case  of  (a)  dry  sand  and  (b)  satu¬ 
rated  sand  created  after  detonation  of  a  100  g 
C4  HE  shallow-buried  mine,  DOB  =  30  mm,  and 
post-detonation  time  =  210  ms.  For  clarity,  the 
steel  barrel  and  the  lower  portion  of  the  sand 
domain  are  not  shown 

and  (b).  It  should  be  noted  that,  while  determination 
of  the  bubble  height  is  fairly  straightforward,  deter¬ 
mination  of  the  bubble  radius  was  associated  with 
some  uncertainty  both  in  the  experiments  and  in  the 
computational  analysis. 

3.2.2  Saturated  sand 

The  ability  of  the  present  material  model  for  sand  to 
account  for  the  main  observations  obtained  during 
experimental  investigation  of  landmine  detonation 
in  fully  saturated  sand  (average  degree  of  saturation 
/So  ~  1.0)  [33]  is  discussed  next.  A  comparison  of  the 
computational  and  experimental  results  pertaining 
to  the  effect  of  SOD  on  the  temporal  evolution  of 
the  sand-overburden  bubble  height  and  radius  are 
shown  respectively  in  Figs  8(c)  and  (d).  The  agree¬ 
ment  between  the  computational  results  and  their 
experimental  counterparts  displayed  in  Figs  8(c)  and 
(d)  is  somewhat  better  than  in  the  case  of  dry  sand, 
Figs  8(a)  and  (b).  This  finding  was  also  expected  since 
the  sand  used  in  reference  [33]  contained  an  undis¬ 
closed  amount  (believed  to  be  between  10  and  20 
per  cent)  of  clay.  Since  the  CU-ARL  visco-plastic  sand 
model  did  not  take  into  account  the  presence  of  clay 
and  since  the  effect  of  clay  is  generally  found  to  be 
more  pronounced  in  dry  sand  (e.g.  reference  [9]),  bet¬ 
ter  computation/ experiment  agreement  was  expected 
in  the  case  of  saturated  sand. 
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Time,  ms  Time,  ms 

(b)  (d) 

Fig.  8  Temporal  evolutions  of  (a)  and  (c)  sand-bubble  height  and  (b)  and  (d)  sand-bubble  radius 
for  (a)  and  (b)  dry  sand  and  (c)  and  (d)  saturated  sand  at  DOB’s  of  30  mm  and  80  mm 
following  detonation  of  a  100  g  C4  HE  shallow-buried  mine 


3.3  Discussion 

The  results  presented  in  sections  3.1  and  3.2  sug¬ 
gest  that  the  CU-ARL  visco-plastic  sand  model  when 
used  in  conjunction  with  the  appropriate  transient 
non-linear  dynamics  simulations  can  reasonably  well 
account  for  the  magnitude  of  the  dynamic  loads  and 
for  spatial  distribution  and  temporal  evolution  of 
the  sand  deformation  accompanying  detonation  of 
shallow-buried  mines  in  dry  or  saturated  sand. 

This  conclusion  is  further  supported  by  the  fact 
that  the  corresponding  experimental  results  were 
associated  with  substantial  uncertainty  (typically  the 


standard  deviation  was  25-30  per  cent  of  the  mean 
value) .  Also,  there  are  additional  phenomena  that  were 
not  accounted  for  in  the  transient  non-linear  dyna¬ 
mics  analysis  of  the  mine  detonation.  For  instance, 
blast  waves,  water  evaporation,  etc. 

The  extent  of  agreement  between  the  CU-ARL  visco¬ 
plastic  sand  model-based  computational  results  and 
the  experimental  results  can  also  be  considered  rea¬ 
sonably  good  considering  the  fact  that  the  soils  used 
in  the  three  sets  of  experiments  had  varying  amounts 
of  clay  and  other  inorganic  and  organic  matter,  as 
well  as  different  average  particle  size  and  particle  size 
distributions. 
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4  VALIDATION  OF  THE  CU-ARL  VISCO¬ 
PLASTIC  SAND  MODEL  IN  THE  LOW 
DEFORMATION-RATE  REGIME 

In  this  section,  the  CU-ARL  visco-plastic  sand  model 
is  tested  by  carrying  out  a  number  of  computational 
transient  non-linear  dynamics  analyses  under  rela¬ 
tively  low  deformation-rate  conditions  and  by  com¬ 
paring  the  computed  results  with  their  experimental 
counterparts.  All  the  calculations  carried  out  in  this 
section  were  also  done  using  ABAQUS/ Explicit  [30]. 

To  validate  the  present  CU-ARL  visco-plastic  sand 
model  in  the  low  deformation-rate  regime,  the  model 
was  first  used  in  the  simulations  of  simple  uniaxial 
stress  relaxation,  uniaxial  creep,  and  triaxial  com¬ 
pression  tests  in  the  same  way  as  it  was  done  in 
reference  [26] .  Since  these  results  are  quite  similar  to 
the  ones  reported  in  reference  [26] ,  and  since  the  tests 
in  question  are  relatively  simple  and  do  not  fully  chal¬ 
lenge  the  present  material  model,  they  are  omitted 
here.  Instead,  the  case  of  a  pneumatic  tire  rolling  on 
dry  or  fully  saturated  sand-based  terrain  is  considered. 
This  problem  is  of  key  importance  to  the  prediction  of 
off-road  vehicle  performance. 

Vehicle  performance  on  unpaved  surfaces  is  impor¬ 
tant  to  the  military  as  well  as  to  agriculture,  forestry, 
mining,  and  construction  industries.  The  primary 
issues  are  related  to  the  prediction  of  vehicle  mobility 
(e.g.  traction,  drawbar  pull,  etc.)  and  to  the  deforma¬ 
tion/  degradation  of  terrain  due  to  vehicle  passage.  The 
central  problem  associated  with  the  off-road  vehicle 
performance  and  terrain  changes  is  the  interaction 
between  a  deformable  tire  and  deformable  unpaved 
road  (i.e.  sand/ soil).  Early  efforts  in  understanding  the 
tire-sand  interactions  can  be  traced  to  Bekker,  at  the 
University  of  Michigan  and  the  US  Army  Land  Loco¬ 
motion  Laboratory  [34-36].  A  good  overview  of  the 
research  carried  out  in  this  area  from  the  early  1960s  to 
the  late  1980s  can  be  found  in  two  seminal  books,  one 
by  Yong  etal.  [37]  and  the  other  by  Wong  [38] .  Over  the 
last  25  years,  there  have  been  numerous  investigations 
dealing  with  numerical  modelling  and  simulations  of 
tire-sand  interactions.  An  excellent  summary  of  these 
investigations  can  be  found  in  Shoop  [39] . 


4. 1  Tire  model  and  validation 

Before  the  problem  of  tire-sand  interactions  can  be 
considered,  a  physically-realistic  tire  model  needs  to 
be  created  and  validated.  Modern  tires  are  structurally 
quite  complex,  consisting  of  layers  of  belts,  plies,  and 
bead  steel  imbedded  in  rubber.  Materials  are  often 
anisotropic,  and  rubber  compounds  vary  throughout 
the  tire  structure.  Models  developed  for  tire  design  are 
extremely  detailed,  account  for  each  material  within 
the  tire,  and  enable  computational  engineering  ana¬ 
lyses  of  internal  tire  stresses,  wear,  and  vibrational 


response.  However,  since,  in  the  present  work,  only 
deformation  of  the  tire  regions  in  contact  with  sand 
and  the  tire’s  ability  to  roll  over  a  deformable  sur¬ 
face  is  of  concern,  a  simpler  model  can  be  employed. 
This  would  provide  for  better  computational  efficiency 
without  significant  loss  in  essential  physics  of  the 
tire-sand  interactions.  Towards  that  end,  a  simpler 
ribbed-tread  tire  model  (described  below)  of  the  type 
often  used  for  harmonic  vibration  modal  analysis  is 
used  here. 

The  tire  considered  in  the  present  work  (the 
Goodyear  Wrangler  HT  235/75  R15)  is  a  modern  radial 
tire,  which  is  composed  of  numerous  components. 
Each  of  these  components  contributes  to  the  struc¬ 
tural  behaviour  of  the  tire  and  is  considered  either 
individually  or  en  masse  when  creating  deformable 
tire  model.  The  tire  direction  convention  used  here  is 
based  on  the  SAE  [40]  standard  definitions  for  vehicle 
dynamics.  Since  the  tire-sand  model  is  used  only  for 
straight-ahead  pure  rolling  (i.e.  zero  slip  angle)  con¬ 
dition,  the  tire  longitudinal  axis  is  in  the  direction  of 
travel  and  the  lateral  or  transverse  direction  is  per¬ 
pendicular  to  travel  and  in  the  horizontal  plane.  Tire 
deflection  at  a  given  level  of  inflation  pressure  is  a  pri¬ 
mary  measure  of  the  tire  structural  response  to  vertical 
load.  Deflection  is  defined  as  the  difference  between 
the  unloaded  and  the  loaded  section  height  and  is  usu¬ 
ally  normalized  by  the  unloaded  section  height  and 
reported  in  percentage.  A  thorough  discussion  of  tire 
mechanics  can  be  found  in  Clark  [41] . 

The  Goodyear  Wrangler  HT  tire  was  evaluated  in 
reference  [39]  for  its  deformation  characteristics  (i.e. 
the  deflection  and  the  contact  area)  on  a  rigid  sur¬ 
face.  The  tire  was  modelled  in  the  present  work  using 
a  ribbed-tread  tire  model  similar  to  that  used  in  refer¬ 
ence  [39] .  The  side-walls  and  the  carcass  are  composed 
of  a  single  layer  of  four-node  reduced-integration  shell 
elements  (S4R)  with  material  properties  representing 
the  composite  behaviour  through  the  carcass  thick¬ 
ness.  The  tread-cap  is  constructed  using  linear,  hybrid 
continuum  elements  (C3D8H) ,  with  constant  pressure 
(simulating  the  nearly  incompressible  nature  of  rub¬ 
ber).  The  result  is  a  reasonable  approximation  of  both 
the  structural  behaviour  and  the  contact  patch.  The 
tire  cross-section  can  be  found  in  Fig.  43  of  refer¬ 
ence  [39] .  To  take  advantage  of  inherent  symmetry  of 
the  model,  tire  model  was  cut  in  half  using  a  verti¬ 
cal  symmetry  plane  along  the  longitudinal  axis.  The 
wheel-rim  and  the  lug-nuts  are  modelled  as  shell  rigid 
bodies.  One  mass  element  was  placed  at  the  wheel’s 
centre  point  on  the  symmetry  plane  and  assigned  a 
mass  consistent  with  the  portion  of  the  vehicle  weight 
supported  by  the  wheel. 

To  validate  the  tire  model  a  multi-step  compu¬ 
tational  analysis  of  tire  contact  and  rolling  over  a 
rigid-road  surface  was  investigated  first.  The  analysis 
included  tire  inflation,  lowering  to  the  road  surface, 
vertical  loading,  and  rolling  (the  wheel’s  centre  point 
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was  translated  while  allowing  it  to  freely  rotate  about 
its  axis  due  to  friction,  simulating  a  towed-wheel  case) . 
The  tire-road  contact  was  assumed  to  be  frictionless 
during  the  loading  step,  and  friction  was  added  during 
the  rolling  step.  During  the  towing  step  the  tire  was 
accelerated  from  rest  to  1  m/s  at  a  constant  accelera¬ 
tion  of  1  m/s2  for  1  s  and  then  it  was  allowed  to  roll  at 
a  constant  velocity  of  1  m/s  for  two  more  seconds. 

A  comparison  of  the  computed  and  measured  tire 
deflections  when  placed  in  contact  with  a  rigid-road 
surface  and  subjected  to  different  vertical  loads  (from 
0  to  8000  N)  at  three  inflation  pressures  (241  kPa / 
35psi,  179kPa/26psi,  and  138kPa/15  psi)  is  shown 
in  Fig.  9(a).  The  suggested  inflation  pressure  for  the 
tire  is  241  kPa.  Lower  inflation  pressures  are  some¬ 
times  used  to  reduce  sinkage  when  driving  in  off-road 


(a) 


Fig.  9  The  effects  of  vertical  load  and  tire-inflation  pres¬ 
sure  on  the  (a)  deflection  and  (b)  contact-patch 
size  for  a  Goodyear  Wrangler  HT  tire  in  contact 
with  a  hard-road  surface 


conditions  and  for  minimizing  damage  to  unpaved 
travel  surfaces.  Consequently,  the  two  lower  inflation 
pressures  were  also  evaluated.  Performance  at  a  range 
of  inflation  pressures  is  also  of  interest  to  industries 
using  vehicles  with  central  tire  inflation  systems  (i.e. 
military,  forestry,  and  agriculture).  The  results  dis¬ 
played  in  Fig.  9(a)  show  that  the  model  predictions 
follow  the  same  general  trend  as  the  experimental  data 
reported  in  reference  [39]  for  all  tire  pressures,  with 
the  largest  difference  occurring  at  the  lower  inflation 
pressure.  This  can  be  readily  justified  since  undesir¬ 
able  tread  and  sidewall  behaviour  can  occur  when  tires 
are  underinflated,  and  these  pressures  and  behaviours 
are  neither  within  the  design  range  nor  accurately 
accounted  for  in  the  model.  On  the  other  hand,  the 
agreement  between  data  and  models  at  the  standard 
inflation  pressure  is  quite  good. 

A  comparison  between  the  model-predicted  and 
the  measured  tire  contact  area  results  is  displayed  in 
Fig.  9(b).  The  contact  areas  are  based  on  the  perime¬ 
ter  of  the  contact,  without  accounting  for  voids  within 
the  area  due  to  tread  design.  In  general,  the  agreement 
between  the  model  and  measured  data  is  reasonable. 

The  distributions  of  the  contact  stresses  over  the 
contact-patch  area  at  different  vertical  loads  and  infla¬ 
tion  pressures  were  also  determined  in  the  present 
work  and  the  results  compared  with  their  experimen¬ 
tal  counterparts  as  reported  in  reference  [39].  Both 
the  computed  and  the  measured  results  revealed  high 
stress  values  at  the  tire  shoulder  and,  to  a  lesser 
extent,  along  the  tire  centre-line.  The  overall  compu¬ 
tation/experiment  agreement  was  quite  good.  Due  to 
copyright  restrictions,  the  measured  contact  stress  dis¬ 
tribution  results  could  not  be  reproduced  here  and, 
hence,  the  corresponding  computed  results  are  also 
omitted. 

As  far  as  the  tire-rolling  step  computational  results 
are  concerned,  they  were  found  to  predict  lower  val¬ 
ues  of  hard- surface  rolling-resistance  when  compared 
to  the  corresponding  experimental  data.  For  example, 
at  a  tire-road  friction  coefficient  of  0.825  (consistent 
with  the  case  of  an  asphalt  pavement),  an  inflation 
pressure  of  241  kPa,  and  a  tire  velocity  of  8  km/h, 
the  computed  rolling-resistance  was  ~10N  while  the 
corresponding  experimental  values  were  ~25  N.  This 
finding  is  consistent  with  the  fact  that  hard-surface 
rolling-resistance  in  tires  contains  a  substantial  contri¬ 
bution  arising  from  the  energy  dissipation  associated 
with  viscoelastic  behaviour  of  rubber  compounds  in 
the  tire,  ffowever,  to  improve  computational  effi¬ 
ciency,  viscoelastic  properties  of  the  tire  materials 
were  not  considered  in  the  present  tire  model.  Conse¬ 
quently,  the  computed  hard-surface  rolling-resistance 
forces  reflect  only  the  contribution  of  the  tire-road 
interfacial  shear.  In  the  case  of  hard- surface  rolling- 
resistance,  the  contributions  of  the  viscoelastic  mate¬ 
rial  behaviour  and  the  interfacial  friction  are  generally 
quite  comparable.  In  sharp  contrast,  when  a  tire 
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rolls  over  deformable  terrain,  the  tire  deforms  signifi¬ 
cantly  less  and,  consequently,  the  rolling-resistance 
is  dominated  by  interfacial  shear.  Thus,  despite  the 
aforementioned  computation/ experiment  discrepan¬ 
cies  regarding  the  hard-surface  rolling-resistance,  the 
present  tire  model  is  deemed  suitable  for  the  analysis 
of  a  tire  rolling  over  a  sand  road. 

Based  on  the  results  presented  and  discussed  in  this 
section,  it  was  concluded  that  the  ribbed-tread  tire 
model  for  the  Goodyear  Wrangler  HT  tire  is  a  good 
compromise  between  physical  reality  and  computa¬ 
tional  efficiency. 

4.2  Analysis  of  tire  rolling  in  sand 

Before  describing  the  details  of  the  finite-element 
analysis  used  to  model  tire  rolling  in  sand,  a  brief 
discussion  is  provided  regarding  tire-sand  interfa¬ 
cial  friction  that  plays  a  major  role  in  the  transfer  of 
forces  parallel  to  the  tire-sand  interface.  Friction  is  a 
complex  phenomenon  that  is  controlled  by  various 
mechanical,  thermal,  physical,  and  chemical  effects 
and,  hence,  is  often  modelled  using  empirical  and 
semi-empirical  laws  and  parameterizations.  These 
laws/parameterizations  generally  apply  only  to  the 
systems  and  conditions  under  which  they  were  deter¬ 
mined  and,  hence,  the  most  reliable  way  to  quantify 
friction  in  a  given  system  is  to  measure  it. 

Typically,  tire-road  interfacial  friction  is  measured 
using  a  vehicle  traction  or  braking  test  and  the 
results  reported  either  as  an  average  value  of  the 
traction/friction  coefficient  or  as  a  traction  force 
versus  wheel-slip  (relative  tire-road  velocity)  curve. 
Both  of  these  could  be  readily  implemented  in 
ABAQUS/ Explicit.  Unfortunately,  the  tire-road  trac¬ 
tion/friction  data  were  not  available  for  the  present 
work.  To  overcome  this  problem,  a  simple  Coulomb 
friction  model  based  on  a  constant  value  of  the  friction 
coefficient  was  used  and  the  tire-sand  analysis  was 
focused  on  the  computation  of  the  rolling-resistance 
forces  (dominated  by  the  deformation  of  sand)  rather 
than  on  the  traction  effort  forces  (dominated  by  the 
interfacial  shear/friction  phenomena). 

4.2. 1  Finite-element  model  of  the  tire  and 
sand-based  terrain 

To  assess  the  ability  of  the  CU-ARL  visco-plastic  sand 
model  to  account  for  the  deformation  behaviour  of 
sand  under  lower  rate  yet  more  complex  loading 
conditions,  a  transient  non-linear  dynamics-based 
computational  analysis  of  the  interaction  between  a 
Goodyear  Wrangler  HT  tire  and  a  deformable  sand- 
based  terrain  was  carried  out  in  the  present  work. 
Since  the  experimental  portion  of  this  work  is  pending, 
a  comparison  between  the  computed  results  and  their 
experimental  counterparts  could  not  be  carried  out. 
Instead,  only  the  computed  results  are  presented 


and  discussed  with  respect  to  the  general  field-test 
observations  regarding  the  tire-sand  interactions. 

The  finite -element  mesh  model  used  in  the 
present  work  to  analyse  the  interactions  between 
the  Goodyear  Wrangler  HT  tire  and  sand-based  ter¬ 
rain  is  displayed  in  Fig.  10.  The  tire  model  used  was 
described  in  section  4.1  and,  hence,  will  not  be  dis¬ 
cussed  here.  As  far  as  the  sand-based  terrain  model  is 
concerned  it  was  represented  using  a  2700  x  1700  x 
450  mm  cuboidal  domain.  The  domain  was  meshed 
using  74  000  eight-node  C3D8R  solid  elements  and 
9200  eight-node  CIN3D8  infinite  elements.  The  later 
elements  were  used  to  model  far-field  sand  regions 
and  avoid  stress-wave  reflection  at  the  sand-domain 
lateral  and  bottom  surfaces.  Finer  30  mm  edge-length 
cubic  C3D8R  elements  were  used  in  the  portion  of  the 
sand  domain  that  was  traversed  by  the  rolling  tire. 
The  sand  domain  containing  C3D8R  elements  was 
filled  with  CU-ARL  visco -plastic  sand  material  while 
the  infinite  elements  were  filled  with  an  ‘elastic’  sand 
material  with  a  Young’s  modulus  and  a  Poisson’s  ratio 
matching  those  of  the  CU-ARL  visco-plastic  sand.  A 
constant  Coulomb  friction  coefficient  was  applied  at 
the  tire-sand  interface. 

Tire-rolling  analysis  was  conducted  by  first  inflat¬ 
ing  the  tire,  then  allowing  the  wheel  to  sink  into  the 
sand  due  to  gravity,  next  accelerating  the  wheel  trans- 
lationally  to  the  desired  speed,  and  then  keeping  it  at 
a  constant  velocity.  The  last  two  steps  were  accom¬ 
plished  by  prescribing  a  longitudinal  velocity  to  the 
wheel  centre  node  which  creates  the  ‘towed-wheel’ 
conditions  (i.e.  the  rotation  of  the  wheel  is  caused  by 
the  tire-sand  interfacial  longitudinal  shear  forces). 

An  example  of  the  results  obtained  in  this  portion 
of  the  work  is  displayed  in  Figs  11  to  13.  The  results 
were  obtained  under  the  following  conditions:  upon 
inflating  the  tire  and  allowing  it  to  sink  into  and  settle 
with  the  sand  below,  the  tire  is  accelerated  longitudi¬ 
nally  at  a  constant  acceleration  of  2777.8  mm/s2  for  1  s 
and  then  its  velocity  was  maintained  at  2777.8  mm/s 
(10  km/h)  for  another  2  s.  Two  cases,  one  involving  dry 


Fig.  10  Computational  domains  used  in  the  rolling 
analysis  of  a  Goodyear  Wrangler  HT  tire  over 
deformable  sand-based  terrain 
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Fig.  1 1  Temporal  evolution  of:  (a)  vertical  sinkage  of 
the  tire  and  (b)  the  motion  resistance  coeffi¬ 
cient  during  tire  roll  over  dry  and  saturated 
sand-based  terrain 


sand  and  the  other  involving  saturated  sand-based 
terrain  were  considered. 

The  results  displayed  in  Figs  11(a)  and  (b)  show 
the  variation  of  wheel  sinkage  (i.e.  the  downward  dis¬ 
placement  of  the  wheeled  centre)  and  the  required 
drawbar  pull  (the  reaction  force  acting  on  the  wheel 
centre  in  the  longitudinal  direction)  with  respect  to 
the  tire-roll  time.  The  results  displayed  in  these  figures 
pertain  to  the  constant-velocity  portion  of  the  simu¬ 
lation.  Furthermore,  the  drawbar  pull  force  displayed 
in  Fig.  1 1  (b)  is  normalized  by  the  normal  force  and 
the  resulting  quantity  termed  the  motion-resistance 
coefficient.  The  results  are  reasonable  since  they  show 
that  a  higher  wheel  sinkage  is  obtained  in  the  case 


Fig.  12  Tire  tracks  left  in:  (a)  dry  and  (b)  saturated 
sand-based  terrain 

of  saturated  sand  (mainly  due  to  lower  associated 
shear  strength)  and  a  higher  value  of  the  drawbar  pull 
(mainly  caused  by  a  larger  volume  of  the  sand  that 
needs  to  be  displaced  by  the  advancing  wheel).  As 
stated  earlier,  due  to  the  use  of  a  simple  Coulomb 
friction  model  at  the  tire-sand  interface  and  a  rela¬ 
tively  small  value  of  the  friction  coefficient  (0.4)  the 
required  draw  bar  pull  is  dominated  by  the  motion 
resistance  (i.e.  by  the  deformation  resistance  of  the 
sand) .  As  mentioned  earlier,  rolling/ motion  resistance 
is  the  sum  of  the  forces  resisting  vehicle  motion.  These 
forces  are  generally  caused  by  the  internal  friction  of 
the  moving  parts  of  the  vehicle,  the  air  drag,  the  inter¬ 
nal  resistance  of  the  tires  due  to  flexing  of  the  belts  and 
plies  and  viscous  properties  of  the  rubber  compounds, 
and  the  added  resistance  due  to  the  deformation  of 
the  road/terrain.  The  internal  resistance  of  the  tire 
was  determined  earlier  by  rolling  it  on  a  hard  surface. 
These  values  were  subtracted  from  the  total  motion 
resistance  values  in  order  to  reveal  the  portion  of  the 
motion  resistance  arising  from  the  deformation  of  the 
sand-based  terrain  and  plotted  in  Fig.  11(b). 

An  example  of  the  track  marks  left  in  the  sand  by  the 
rolling  tire  is  displayed  in  Figs  12(a)  and  (b)  for  dry  and 
saturated  sand-based  terrains,  respectively.  Clearly  the 
track  depth  is  larger  in  the  case  of  saturated  sand  and 
the  track-side  bulldozing  effect  is  more  pronounced. 
These  observations  are  consistent  with  typical  field- 
test  observations  (e.g.  reference  [42]). 

Contour  plots  showing  the  distribution  of  the  maxi¬ 
mum  shear  stress  in  the  sand  beneath  the  tire  is 
displayed  in  Figs  13(a)  and  (b).  In  both  the  case 
of  dry  sand  and  saturated-sand  based  terrain,  the 
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Fig.  13  Typical  spatial  distribution  of  the  maximum 
shear  stress  in  sand  beneath  the  advancing 
wheel:  (a)  dry  and  (b)  saturated  sand-based 
terrain 


largest  values  of  the  maximum  shear  stress  are  point¬ 
ing  downwards  and  away  from  the  wheel,  revealing 
the  location  where  a  new  shearing  surface  begins  to 
develop.  These  observations  are  quite  consistent  with 
experimental  observations  (e.g.  reference  [42]). 

Based  on  the  results  presented  in  this  section,  it 
can  be  concluded  that  qualitatively  the  CU-ARL  visco¬ 
plastic  sand  model  can  quite  well  account  for  typical 
observations  made  in  field-tests  of  the  off-road  vehi¬ 
cle  travel.  Considering  the  fact  that  the  sand  model 
was  also  validated  quantitatively  relative  to  the  simple 
laboratory  uniaxial  creep  and  stress  relaxation  tests, 
and  quasi-static  three-dimensional  compression  tests, 
it  can  be  concluded  that  the  model  is  quite  promising. 
In  our  future  communications,  a  more  detailed  quan¬ 
titative  comparison  will  be  given  between  the  com¬ 
putational  results  based  on  the  current  sand  model 
and  their  experimental  counterparts  relative  to  the 
tire-sand  interactions. 


5  SUMMARY  AND  CONCLUSIONS 

Based  on  the  results  obtained  in  the  present  work, 
the  following  main  summary  remarks  and  conclusions 
can  be  drawn. 

1 .  The  visco-plastic  material  model  recently  proposed 
by  Tong  and  Tuan  [26J  has  been  expanded  to  take 
into  account  the  effect  of  moisture  and  the  effect  of 
rate  dependency  of  the  bulk  modulus. 


2.  The  new  model,  named  the  CU-ARL  visco-plastic 
sand  model  has  been  tested  against  a  series  of 
laboratory  experiments  and  field-test  data  and 
observations,  both  in  the  high-deformation-rate 
(e.g.  detonation  of  mine  shallow-buried  in  sand 
and  soil-ejecta-target-structure  interactions)  and 
in  the  low-deformation-rate  regime  (e.g.  sinkage 
and  rolling  of  a  pneumatic  tire  over  deformable 
sand-based  terrain). 

3.  In  the  high-deformation-rate  regime,  the  present 
sand  model  has  been  found  to  clearly  outperform 
the  Laine  and  Sandvik  [5]  sand  compaction  model. 
The  latter  model  is  very  frequently  used  in  the 
computational  blast  analysis  community. 

4.  In  the  low  deformation-rate  regime,  the  CU-ARL 
visco -plastic  sand  model  has  been  found  to  yield 
predictions  regarding  the  extent  of  wheel  sinkage, 
drawbar  pull,  and  the  distribution  of  stresses  in 
sand  beneath  the  wheel  which  are  fully  consistent 
with  the  corresponding  field-test  observations  and 
measurements. 
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